Multilocus Data Analysis Reveal the Diversity of Cryptic Species in the Tillandsia ionantha (Bromeliaceae: Tillansiodeae) Complex

Independent evolutionary lineages or species that lack phenotypic variation as an operative criterion for their delimitation are known as cryptic species. However, these have been delimited using other data sources and analysis. The aims of this study are: (1) to evaluate the divergence of the populations of the T. ionantha complex; and (2) to delimit the species using multilocus data, phylogenetic analysis and the coalescent model. Phylogenetic analyses, genetic diversity and population structure, and isolation by distance analysis were performed. A multispecies coalescent analysis to delimit the species was conducted. Phylogenetic analysis showed that T. ionantha is polyphyletic composed of eight evolutionary lineages. Haplotype distribution and genetic differentiation analysis detected strong population structure and high values of genetic differentiation among populations. The positive correlation between genetic differences with geographic distance indicate that the populations are evolving under the model of isolation by distance. The coalescent multispecies analysis performed with starBEAST supports the recognition of eight lineages as different species. Only three out of the eight species have morphological characters good enough to recognize them as different species, while five of them are cryptic species. Tillandsia scaposa and T. vanhyningii are corroborated as independent lineages, and T. ionantha var. stricta changed status to the species level.


Introduction
Cryptic species refers to taxa that are erroneously lumped into a single nominal species due to the paucity of conspicuous morphological differences, but have been delimited from other sources of data and analysis [1][2][3]. Historically, morphological criteria have been used to identify and delimit species. However, relying solely on morphology will likely result in an underestimation of the number of species [4]. This fact has required the use of other data sources (e.g., molecular, ecological, and biogeographical) and analytical methods (e.g., phylogeographical, population genetics, phylogenetic, and coalescent based) for the recognition of cryptic lineages [2,[5][6][7].
The emergence of molecular and genomic datasets, as well as contemporary species concepts such as the one that considers species as independent evolutionary lineages [8,9], have brought species delimitation to an interesting crossroads with various methodological and philosophical approaches [10]. Currently, several effective techniques and methods have been developed for the more accurate delimitation and classification of species; among these are phylogenetic and non-phylogenetic methods [11] as well as Bayesian methods based on the coalescent model [12][13][14]. Recently, the application of the coalescent theory provides a fundamentally different and stronger framework to objectively identify cryptic and/or allopatric species using genetic data [10]. In addition, coalescent methods allow us to infer dynamics of divergence, interaction of evolutionary processes and relationships among taxa [15].
In this study, our model taxa are the so-called Tillandsia ionantha Planch. complex. Ancona et al. [16] postulated that the T. ionantha complex is composed of two species, five varieties, and two forms: Tillandsia scaposa (L.B. Sm.) Ehlers  The complex is widely distributed in the tropical dry forests, savannahs, and oak forests in transition with the dry forests of Mexico and Central America, and its elevation ranges from 60 to 1600 m. Morphologically, the members of the T. ionantha complex are characterized by having relatively small rossetes compared to others members of the subgenus Tillandsia (20 cm length). In general, they are epiphytic and occasionally rupicolous. They present 10 (-20) cm tall rosettes, acaulescent or caulescent, with green or greyish leaves that turn red at anthesis and that are covered by a silverylepidote coating; they have 2-3 (4) flowers nested in the center, with violet petals bearing exerted stamens and stigmas [17].
Tillandsia scaposa, was originally described as a variety of T. ionantha (=T. ionantha var. scaposa L.B.Sm.). Ehlers [18] nowadays treated as a different species due to the presence of a pedunculated inflorescence. Tillandsia ionantha var. vanhyningii was described by Foster [19]. At first, this species was considered as a new species, but when compared with T. ionantha, no differences were found other than the rupicolous habit and caulescent rosettes, but the arrangement and position of the flowers, color of the petals and exserted stamens and pistil led him to conclude that it was simply a variant. After the publication of Ancona et al. [16], Beutelspacher and García-Martínez [20] changed the status of T. ionantha var. vanhyningii to species. Other varieties and forms described, do not present morphological variation in reproductive structures; the descriptions were made based on vegetative morphological characters, which are highly variable between and within wild populations [16,21]. However, among these taxa there are populations that are isolated by geographical barriers that might play an important role in allopatric speciation.
The phylogenetic position of the Tillandsia ionantha complex has been evaluated using molecular data [22]. With four cpDNA markers, the results showed that T. ionantha in its current constituency is polyphyletic. Two large clades were also detected: one integrated by the populations of the Pacific slope and another by the slopes of the Gulf of Mexico and Central America; in addition to T. scaposa and T. ionantha var. vanhyningii form, two clades isolated from the rest of the varieties of T. ionantha [22]. It is important to highlight that Rodríguez's objective was to determine the phylogenetic position of the T. ionantha complex within the subgenus Tillandsia and not to delimit species in the complex. But his results suggest lineage hypotheses that should be tested. Thus, in this study, we are interested in answering the following questions: Are Tillandsia scaposa and T. vanhyningii independent species? Is there population genetic structure and genetic isolation among populations? If such a divergence exists, is geographic distance a factor in the evolution of T. ionantha populations? Are cryptic species present in the T. ionantha complex? Therefore, the aims of this study are (1) to evaluate the genetic variation between the populations and varieties of the T. ionantha complex, and (2) to delimit its species using multiple genomic regions and phylogenetic analyses.

Phylogenetic Analysis
The topology of the phylogenetic analysis with Bayesian inference and Maximum Likelihood (ML) were identical in both sets of markers. Here we only present the phylogenetic trees of Bayesian inference for cpDNA and nDNA markers. On the branches, we feature the bootstrap support values resulting from the phylogenetic analysis based upon ML. The phylogenetic tree obtained from the analysis of three cpDNA markers showed that Tillandsia ionantha in its current circumscription is polyphyletic (Figure 1). In general, the complex is arranged in three major clades: (1) a clade that includes the T. ionantha lineage 1, T. ionantha var. stricta, T. scaposa, T. ionantha lineage 2 plus the species T. strepotphylla Scheidw. ex E. Morren, T. pruinosa Sw. and T. seleriana Mez. (2) An intermediate clade that includes the lineage 3 and lineage 4. These populations are phylogenetically isolated from the rest of the varieties of T. ionanatha and more closely related to the species of the outgroup from T. lydiae Ehlers to T. variabilis Schltdl. (3) Finally, a clade that includes T. ionantha lineage 5 and T. dasyliriifolia Baker as a sister species. Other outgroup species plus T. vanhyningii are inserted into this clade.
Within the complex, eight evolutionary lineages with high support values are observed. The populations of Tillandsia ionantha var. stricta, T. vanhyningii, and T. scaposa are recovered as three monophyletic lineages, whereas populations of T. ionantha var. ionantha segregated into five lineages. Lineage 1 is made up of the Central American populations and includes the type locality of T. ionantha var. zebrina (Nicaragua, Guastatoya and Zacapa). Lineage 2 is made up of the populations of T. ionantha var. ionantha from northeastern Mexico (Jalcomulco, El Higo, Tolimán, and El Troncón). Lineages 3 and 4 comprise the populations of T. ionantha var. ionantha from Pochotitan and Matatán, respectively. Lineage 5 is composed of the populations of T. ionantha var. ionantha from the Pacific lowlands and Lower Balsas Basin, including the type locality of T. ionantha var. maxima (Jalisco, Pilcaya, Santo Tomás de los Plátanos, el Puente, and Huamelula).
Phylogenetic analysis with Bayesian Inference using the PHYC marker array shows a different typology than that obtained with cpDNA markers. With this nuclear marker, T. ionantha in its current district is paraphyleticand sister to T. scaposa. However, within the complex, the eight lineages recognized with the cpDNA markers are recovered with high values of posterior probability and bootstrap ( Figure 2). Tillandsia ionantha lineage 5 and T. vanhyniningii are sister lineages.

Genetic Diversity and Structure of cpDNA
The alignment was 3367 bp long and consisted of 229 nDNA sequences, with 45 polymorphic sites resulting in 13 haplotypes. The estimated haplotype diversity was Hd = 0.83 and the nucleotide diversity π = 0.0025. Genetic diversity across all populations (hT = 0.888; vT = 0.884) was greater than within population average value (hS = 0.076; vS = 0.009). The 13 haplotypes identified are strongly associated with biogeographic districts representing the same pattern observed with phylogenetic analyses ( Figure 3A,B). The Oaxacan Plateu biogeographic district featured the highest number of haplotypes (haplotypes 8, 9, 10, and 11), followed by the Sinaloan district with three haplotypes (haplotypes 5, 6 and 7). On the other hand, haplotype 4 is the one with the widest geographical distribution and it is located in the Nayarit-Guerrero, Lower Balsas Basin, and Tehuanan biogeographic districts. Haplotype 13 has the second widest geographical distribution, being found in two Central American biogeographic districts, El Mosquito and Tapachultecan. The rest of the haplotypes are located in a single biogeographic district ( Figure 3A,B).  The presence of private haplotypes reflects a high genetic differentiation (F ST = 0.98) within cpDNA. On the other hand, the differentiation value considering distances between haplotypes (ordered alleles) was N ST = 0.989 and significantly different from the G ST value (0.914; p < 0.05), which denotes a strong phylogeographic structure. Analysis of molecular variance (Table 1) showed that 99% of the variation is located between populations and 1% of the variation is located within populations. When evaluating the genetic differentiation between the six varieties of Tillandsia ionantha, a negative, non-significant value of F CT was estimated, meaning that there is no variation between the populations of the varieties within the current circumscription of T. ionantha. The variation among the eight lineages is the  (Table S1). In Figure 4A the strong genetic structure between the populations is shown. There are populations that behave as a single genetic population, for example the NIC, GUAS and ZAC populations, are strongly differentiated from the rest of the populations. The Mantel test showed that the genetic differentiation of cpDNA between populations is correlated with geographic distance (R2 = 0.272, p = 0.0001).

Genetic Diversity and Structure of cpDNA
The alignment was 3367 bp long and consisted of 229 nDNA sequences, with 45 polymorphic sites resulting in 13 haplotypes. The estimated haplotype diversity was Hd = 0.83 and the nucleotide diversity π = 0.0025. Genetic diversity across all populations (hT =

Genetic Diversity and Structure of nDNA
The alignment was 1114 bp long and consisted of 227 nDNA sequences, with 55 polymorphic sites resulting in 88 haplotypes. The estimated haplotype diversity was Hd = 0.95 and the nucleotide diversity π = 0.0047. Genetic diversity across all populations (hT = 0.972; vT = 0.997) was greater than within population average value (hS = 0.737; vS = 0.244). The 88 haplotypes identified grouped into eight phylogroups strongly associated with biogeographic districts, representing the same pattern observed with cpDNA (Figure

Genetic Diversity and Structure of nDNA
The alignment was 1114 bp long and consisted of 227 nDNA sequences, with 55 polymorphic sites resulting in 88 haplotypes. The estimated haplotype diversity was Hd = 0.95 and the nucleotide diversity π = 0.0047. Genetic diversity across all populations (hT = 0.972; vT = 0.997) was greater than within population average value (hS = 0.737; vS = 0.244). The 88 haplotypes identified grouped into eight phylogroups strongly associated with biogeographic districts, representing the same pattern observed with cpDNA ( Figure 3A,C). The biogeographic districts that present a greater diversity of haplotypes were Oaxaca Plateau, Nayarit-Gerrero + Lower Balsas Basin + Tehuanan, Deciduous Forest of Northern of Veracruz and Tapachultecan + Mosquito. The distribution of these haplotypes is star-shaped; which means, low-frequency haplotypes are derived from higher-frequency haplotypes, separated by a single mutational step.
The presence of private haplotypes reflects a high genetic differentiation (F ST = 0.75). The differentiation value considering distances between haplotypes (ordered alleles) was N ST = 0.755 and significantly different from the G ST value (0.242; p < 0.05), which denotes a strong phylogeographic structure. Analysis of molecular variance (Table 1) was similar to that estimated with cpDNA. When analyzing the populations as a single group, it resulted in that 76.59% of the variation is between the populations and only 23.43% of the variation is located within the populations. When evaluating the genetic differentiation between the six varieties of T. ionantha, a non-significant low value of F CT = 0.38 was observed. That is, there is no variation among the populations of the varieties of the current circumscription of T. ionantha using the PHYC marker. Like the cpDNA, the variation among the eight lineages presented the highest F CT = 0.78 value. Pairwise F ST values ranged from 0.0 to 0.936 (Table S2). However, populations that behave as a single genetic unit and strongly differentiated from the rest of the populations were also observed ( Figure 4B). The Mantel test showed that the genetic differentiation of nDNA between populations is correlated with geographic distance (R2 = 0.283, p = 0.010).

Multispecies Coalescent
The multispecies coalescent analysis supported ten distinct clades (NCclusters = 10): two correspond to the outgroup species and the remaining eight are within the Tillandsia ionantha complex and are consistent with the phylogenetic analysis. Figure 5 shows the tree of minimal clusters with high PP values in each of the detected lineages. The gene trees recover as highly supported the eight clades identified with the phylogenetic analyzes ( Figures S1-S4). However, support between lineages is moderate.
test showed that the genetic differentiation of nDNA between populations is correlated with geographic distance (R2 = 0.283, p = 0.010).

Multispecies Coalescent
The multispecies coalescent analysis supported ten distinct clades (NCclusters = 10): two correspond to the outgroup species and the remaining eight are within the Tillandsia ionantha complex and are consistent with the phylogenetic analysis. Figure 5 shows the tree of minimal clusters with high PP values in each of the detected lineages. The gene trees recover as highly supported the eight clades identified with the phylogenetic analyzes ( Figures S1-S4). However, support between lineages is moderate.

Discussion
In the Tillandsia ionantha complex, gene trees using three cpDNA regions and a nuclear marker evidenced distinct evolutionary histories in terms of their phylogenetic position within the subgenus but not at the population level, where populations remain reciprocally monophyletic. Several studies have suggested that topological inconsistency between cpDNA and nDNA markers might due to various factors, including convergent evolution, introgression after hybridization, incomplete lineage sorting, and horizontal transfer [23][24][25][26]. For our purposes, and aiming to delimit putative species, both sets of markers identify the same eight lineages with high support values within the T. ionantha complex. Three of these eight evolutionary lineages, the T. ionantha var. stricta, T. scaposa and T. vanhyningii maintain their morphological, genetic and biogeographic identity. The remaining lineages maintain their genetic and biogeographical identity, but are morphologically cryptic. These results correspond to the allopatric speciation model proposed by Harrison [27], where the process begins with the interruption of gene flow between populations, and then different alleles (or haplotypes) are fixed in the populations, which may or may not acquire different phenotypic characters until they become reciprocally monophyletic. Therefore, these lineages must be considered as evolutionarily independent lineages [8,9].
The diversity and distribution of the haplotypes, the analysis of molecular variation, and the F ST values for the pairwise comparison among populations with both sets of markers are consistent with the eight clades (lineages) obtained with the phylogenetic analysis. The Mantel test on both sets of markers correlates positively with the geographic distance of the populations, supporting the hypothesis of the allopatric speciation related with the geographical barriers and orographic features of Mexico and Central America that are quite diverse [28][29][30], as these factors lead to a great diversity in climates, vegetation types, and microenvironments [31,32] that allow the geographic isolation of populations and therefore genetic isolation, divergence and the speciation of populations (ecological/allopatric Plants 2022, 11, 1706 9 of 17 speciation). For Quercus [33], Agave [34], Bursera [35], Bakerantha [36], and the Tillandsia utriculata complex [37] diversification and divergence of species was suggested to be caused by isolation, expansion, and colonization of new environments.
With our results, we can say that genetic differentiation and population divergence in the T. ionantha complex are affected by geographic distance and geographic barriers. Despite the geographical distance, no morphological changes were observed in the reproductive structures [21], which could suggest a niche conservatism in the populations. In contrast to the vegetative traits, species of genus Tillandsia have developed a great diversity of floral traits that attract a wide variety of pollinators, including insects, birds, and bats [38][39][40]. However, floral morphology in the T. ionantha complex lineages is identical and flowering times overlap; which indicates that they must share the same pollination syndromes and their pollinators must be local and delimited by the same geographical barriers. If so, this hypothesis only supports that geographic barriers prevent the arrival and spread of new individuals/seeds (cpDNA) and pollen (nDNA) from the population of one lineage to the population of another lineage, thus conserving the local gene pool of each lineage. In addition, the high genetic differentiation (F ST ) between the lineages of the T. ionantha complex reflect the high rates of gene flow between the populations of the same lineage, thus promoting rapid monophyly at the species level, indicating that gene flow is beneficial for the maintenance of these lineages and their cohesive evolution [41]. In contrast, if geographic barriers were not operative as isolation factors, gene flow would increase between lineages, preventing haplotype fixation and hampering genetic differentiation between detected lineages [41,42].
Multispecies coalescent methods for species delimitation have been highlighted as more objective, robust, and repeatable than traditional morphological taxonomy. However, the use of such methods for the delimitation of species and, later, to make taxonomic decisions is still widely discussed and hesitantly accepted among botanists, especially when cryptic species are discovered [43]. The multispecies coalescent analysis applied in the present work is consistent with the phylogenetic analyses, population structure, and genetic isolation to delimit eight lineages. Three of these eight lineages have morphological characters that, in practical terms, are easy to recognize in the field, while there is no morphological differentiation so far for the remaining five lineages. The cryptic speciation has been increasingly observed in neotropical groups of angiosperm plants [43][44][45], fungi [46][47][48] and animals [49][50][51], suggesting that biodiversity might be underestimated in this region.
Regarding Tillandsia, this work is the first to using multilocus data, phylogenetic and phylogeographic analyzes, genetic divergence and differentiation, and the coalescent model to delimit species. It is also perhaps the first work that reports the presence of cryptic species in the genus Tillandsia. More work at this same population level will help to identify the presence of more cryptic species in Tillandsia. In Bromeliaceae, other studies [52] used multilocus data and phylogenetic and phylogeographic analysis to delimit the two species in the Tillandsia capillaris Ruiz and Pav. complex. Romero-Soler et al. [36] used multilocus data, phylogenetic analysis, and the coalescent model to delineate species in the genus Bakerantha L.B. Sm. The only two studies that report the presence of other cryptic species in Bromeliaceae are Santos-Leal et al. [53], who through multilocus data, phylogeographic analysis, and coalescent models to identify the presence of a cryptic species in the South American species Pitcairnia lanuginose Ruiz and Pav. Goncalves-Oliveira et al. [54] using microsatellite data and phylogeographic analysis reported several cryptic species in Encholirium spectabile Mart. ex Schult. f.
Under the general concept of lineage species, which differentiates species as "lineages of metapopulations that evolve separately" [9], in this study we regard the results of the data based on trees, morphological discontinuity [21] and geographic and ecological distribution as providing evidence strong enough to recognize the evolutionary independence between different species of the Tillandsia ionantha complex. The current circumscription of the populations and varieties of T. ionantha should change according to phylogenetic, genetic isolation and biogeographic evidence here demonstrated. The morphological variation is evident for three lineages within the T. ionantha complex. The first is T. scaposa that had already been recognized morphologically and its status changed by Ehlers [18] due to the presence of a minute and pedunculated inflorescence, in addition to being ecologically different as it inhabits pine-oak forests in the Guatemalan highlands. Tillandsia vanhyninigii differs morphologically from the other varieties due to the caulescent rosette that can measure up to 20 cm in length, the sheath almost the same size as the blade; furthermore, it inhabits rocky walls at the Sumidero Canyon, Chiapas. The populations of the T. ionantha var. stricta clade are distinctive within the complex by its rosettes, which are compressed and the smallest within the complex, measuring 5-8 (-10) cm [21]; it is restricted in distribution to the Oaxacan Plateau Biogeographic district.
On the other hand there are five lineages that are morphologically indistinguishable, and that we propose to recognize as cryptic species. Three of these lineages had been recognized as varieties of T. ionantha. However, the diagnostic characters employed by some authors are variable within and between populations of the different lineages [21]. The diagnostic character of T. ionantha var. maxima (included in lineage 5) is its large rosettes, which can be twice as large as any other variety of T. ionantha. But the size of the rosette is not taxonomically informative in differentiating the evolutionary lineages found in our analysis. The size of the rosettes in other varieties varies between 5 and 12 cm within the same population. If a 5 cm rosette is considered, the double would be 10 cm and this value still falls within the taxonomic varieties of T. ionantha and does not exactly diagnose T. ionantha var. maxima. The diagnostic character in T. ionantha var. zebrina (included in lineage 1) is the banded pattern of the leaves. This character may or may not be present in individuals within the same population. In addition, we have observed this character in individuals of the T. ionantha lineage 2 from Northeastern Mexico, and therefore this character is not taxonomically informative to delimit these lineages. Finally, there are two lineages in the Sinaloan biogeographic district, being morphologically similar to the rest of the cryptic species within this complex.

Sampling Strategy
In order to know the distribution area of the varieties of the T. ionantha complex and to select the sampling populations, the type specimens and other specimens deposited in the herbaria CICY, K, F, GH, GOET, M, MEXU, MO, NY, P, SEL, US and WU were revised [55]. The revision of the herbaria allowed us to recognize that the varieties T. ionanatha var. maxima and T. ionantha var. zebrina are only known from the type collections (holotypes, isotypes and paratypes). For T. vanhyningii. T. ionantha var. stricta, and T. scaposa we found several specimens but they are restricted to biogeographical areas ( Table 2); while T. ionantha var. ionantha has the widest geographical distribution ( Figure 6). Nineteen populations, each with 10 to 21 individuals were selected (Table 2), leaves were collected and stored in silica gel. The selected populations included the type localities for each variety of T. ionantha, T. vanhyningii and T. scaposa. Similarly, we selected populations to cover the whole geographical distribution of T. ionantha var. ionantha (Figure 6).

DNA Isolation, Amplification, and Sequencing
DNeasy Plant Mini Kit (QIAGEN) was used to extract DNA from 240 individuals from the 19 populations (Table 2) and 35 outgroup species [22], The sequences will be deposited in the Genbank (if you want these sequences contact the corresponding author). The outgroup species selection was based on previous phylogenetic studies of Barfuss et al. [56], and Pinzón et al. [37].

DNA Isolation, Amplification, and Sequencing
DNeasy Plant Mini Kit (QIAGEN) was used to extract DNA from 240 individuals from the 19 populations (Table 2) and 35 outgroup species [22], The sequences will be deposited in the Genbank (if you want these sequences contact the corresponding author). The out-group species selection was based on previous phylogenetic studies of Barfuss et al. [56], and Pinzón et al. [37].

Phylogenetic Analyses
To explore the variability of each marker, we calculated the number of variable and potentially parsimony informative sites using MEGA11: Molecular Evolutionary Genetics Analysis version 11 [59]. Two phylogenetic reconstructions analyses were performed using MrBayes 3.2.7 software [61], one with the cpDNA matrix data and the other with the nDNA matrix data. For each matrix data we used the estimated nucleotide substitution model in MegaX (Table S3). We performed two independent runs with 10,000,000 generations for each matrix data, sampling every 1000 generations. Each run consisted of one cold and three hot chains, unbinding the "tratio", "statefrq" and "shape" parameters on all partitions, as well as setting "rate" as a variable for all partitions. A 50% majority rule consensus tree was constructed after discarding the initial 25% of trees. For the Maximum Likelihood approach, we ran 10 replications to retrieve the best topology in RAxML 8.2.12 [62] and a 1000 replications Bootstrap for nodal support using the CIPRES Science Gateway 3.3 [63], with all other parameters as default. The trees obtained were visualized and edited using FigTree [64].

Genetic Diversity and Structure Analyzes
Like the phylogenetic analyses, we performed the genetic structure analyses separately for both sets of markers. The number of haplotypes (H), the diversity of haplotypes (Hd), the diversity of nucleotides (Π) and polymorphic sites were estimated using the DNSpV.6 program [65]. As each of the eight evolutionary lineages found in the phylogenetic analyzes are restricted to specific biogeographical districts. In this section, we use PopART version 1.7 [66] and the TCS network method [67] to analyze the evolutionary relationships between haplotypes and biogeographic districts where each of the studied populations are distributed (Table). Haplotype distribution patterns were plotted on a map. The indices for ordered (vS, vT) and unordered (hS, hT) haplotypes and differentiation parameters N ST , G ST [67] and N ST [68,69] were estimated using PERMUT v.1.0 [69]. The significant difference between the N ST and G ST values was tested with 10,000 permutations. A significant value of N ST indicates phylogeographic structure and that the most similar haplotypes are geographically close to each other.
We evaluated the population structure with pairwise F ST values calculated with ARLEQUIN version 3.5 [70] and PAST version 4.08 [71] was used to represent the values in a matrix plot. Analysis of molecular variance (AMOVA) [72] were also performed using ARLEQUIN v.3.5 to test the structure between populations using both data sets separately and different hierarchical levels: eight lineages detected by phylogenetic analysis and the current taxonomic classification of the T. ionantha complex [17]. To determine whether divergence among populations is an effect of isolation by distance, we correlated pairwise genetic distance (measured as F ST ) versus geographical distance matrices by a Mantel test using the Software GenAlEx 6.5 [73].

Species Delimitation Using STACEY
A Bayesian coalescent method of species delimitation was performed using STACEY (species tree estimation using DNA sequences from multiple loci) [74], which is an extension of DISSECT [75]. STACEY was implemented in BEAST ver. 2.4.4 [76][77][78].
We ran STACEY using a multilocus dataset (ptDNA and nDNA) with all populations of the T. ionantha complex and two outgroup species. Input xml files were prepared in BEAUti v. 2.4.7 [76] and the corresponding substitution models were set according to results from MegaX (Table S3). The analysis were run for 100 million generations of the MCMC chains, sampling every 5000 generations. Convergence of the stationary distribution was checked by visual inspection of plotted posterior estimates using Tracer ver. 1.7 [79]. After discarding the first 1000 trees as burn-in, the samples were summarised in the maximumclade credibility tree using TreeAnnotator ver. 1.6.1 [77] with a posterior probability limit of 0.5 and summarising of mean node heights. The results were visualised using FigTree ver. 1.3.1 [64].

Conclusions
Multilocus data and different analysis methods support species delimitation and cryptic diversity in T. ionantha. Morphological and genetic evidence corroborates the hypothesis of the evolutionary lineage of Tillandsia scaposa, T. vanhyningii and T. ramireziana (=T. ionantha var. stricta). On the other hand, the remaining five lineages are supported by genetic evidence, but not by morphology. Even though these five cryptic species do not have unambiguous macromorphological characters, the hypothesis that they are populations that have speciated is supported by the occurrence of clear phylogenetic structure (high levels of divergence between the long branches of the phylogenetic trees and high values of support in nodes), geographical and ecological differences. However, for practical reasons and identification of the populations in the field, we will avoid describing and naming them until we have more morphological evidence. Further studies can explore micromorphological and anatomical characters to corroborate these hypotheses of cryptic lineages in T. ionantha.

Taxonomy
Tillandsia scaposa is validated as an evolutionary lineage as proposed by Ehlers [17]; and two varieties of T. ionantha are treated at the species level: T. vanhyningii and T. ramireziana. Likewise, it is proposed that these taxa be eliminated from the T. ionantha complex since the morphological and molecular evidence support this. In such a way that the complex is currently composed of five morphologically cryptic species. Note. Tillandsia ramireziana is a replacement name for T. stricta, since the name already exists to name for species that is distributed in Brazil (Tillandsia stricta Sol. ex Sims). This name is restricted to the populations that are distributed in the Oaxacan Plateau and the Chiapas Highlands (Figure 7). The vegetation is a low deciduous forest type in transition with oak forests in elevation from 200 to 1600 m. Within this species, Koide [80] described two morphological forms. However, our analysis does not recognize the delimitation of both forms as lineages; therefore both names are recognized as synonyms of T. ramireziana.
Eponymy. The specific epithet of this species is in honor to Dra. Ivón Ramírez Morillo for her long career in the taxonomic and systematic studies on Bromeliaceae. Tillandsia

Supplementary Materials:
The following supporting information can be downloaded www.mdpi.com/xxx/s1. Figure S1. Gene tree derived from StarBEAST analysis using the PH marker sequences. Figure S2. Gene tree derived from StarBEAST analysis using the trnT-L-F mar sequences. Figure S3. Gene tree derived from StarBEAST analysis using the rps16-trnQ marker quences. Figure S4. Gene tree derived from StarBEAST analysis using the rpl16-rps3 marker quences. Figure S5. Tree of maximum credibility of the BEAST analysis using three cpDNA mol ular markers (trnT-L-F, rps16-trnQ and rpl16-rps3) and the nDNA PHYC marker and each popu tion as a putative lineage. Table S1. FST values for pairwise comparison between populations of T. ionantha complex using cpDNA sequences. Table. S2. FST values for pairwise comparison tween populations of the T. ionantha complex using nDNA sequences. Table S3. List of sequen generated in this study and deposited in GenBank. Table S4. Names and sequences of the prim used in this study and temperature conditions for PCR. * = Primers used in the amplification; Primers used in sequencing. References [22,57,81,82].  Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/plants11131706/s1. Figure S1. Gene tree derived from StarBEAST analysis using the PHYC marker sequences. Figure S2. Gene tree derived from StarBEAST analysis using the trnT-L-F marker sequences. Figure S3. Gene tree derived from StarBEAST analysis using the rps16-trnQ marker sequences. Figure S4. Gene tree derived from StarBEAST analysis using the rpl16-rps3 marker sequences. Figure S5. Tree of maximum credibility of the BEAST analysis using three cpDNA molecular markers (trnT-L-F, rps16-trnQ and rpl16-rps3) and the nDNA PHYC marker and each population as a putative lineage. Table S1. FST values for pairwise comparison between populations of the T. ionantha complex using cpDNA sequences. Table. S2. FST values for pairwise comparison between populations of the T. ionantha complex using nDNA sequences. Table S3. List of sequences generated in this study and deposited in GenBank. Table S4. Names and sequences of the primers used in this study and temperature conditions for PCR. * = Primers used in the amplification; + = Primers used in sequencing. References [22,57,81,82].